function [FinalCurve]=RibCageApproximation(I,Thres);
%Draw a  closed curve which approximates the cross
%section of the rib cage in CT slices.
%      [FinalCurve]=RibCageApproximation(I,Thres);
%       Input parameters
%       Img=> the CT slice
%       Thres=> the desired threshold for the ribs (bones) extraction.

Img=squeeze(I);

%apply  threshold
BinaryImage=Img>Thres;
%perform morphological open operation
BinaryImage = imopen(BinaryImage,strel([1 1;1 1]));

%find the pixels of the rib cage.
[x,y]=find(BinaryImage==1);

%define the center of the image.
ImgCenter =round( [size(BinaryImage)]./2);


%Estimate the vector between bone pixels and image centre
vector = [x-ImgCenter(1),y-ImgCenter(2)];
%Vector's angle in degrees.
Angle = mod(360+180/pi*atan2(vector(:,2),vector(:,1)),360);


%In each directions (0, 10,...., 360 degrees), define the bone pixel which is closer to the image
%center and keep it for the interpolation process.

Theta=[0:10:360];k=0;

Ind=cell([length(Theta)-1,1]);
for i =1:length(Theta)-1
    [Ind{i,1}] = find(Angle>=Theta(i) & Angle<Theta(i+1));
    if ~isempty(Ind{i,1})
      IndMinDist = FunctionMinDistance(x(Ind{i,1}),y(Ind{i,1}),ImgCenter(1),ImgCenter(2));
      k=k+1;
      IndOfMinDistance(k,1)=Ind{i,1}(IndMinDist(1),1);
       
    end 
end


%Set points with the minimum distance from the center of the ribs cavity in an order
[xi, yi]=PixelsInOrder(x(IndOfMinDistance),y(IndOfMinDistance),ImgCenter(1),ImgCenter(2));

%Cubic interpolation
[xii, yii] = Curvefitting(xi,yi,0.2);
xii = round(xii);yii=round(yii);

%Draw the final curve
FinalCurve=zeros(size(BinaryImage));

for i =1:length(xii)
    FinalCurve(xii(i),yii(i))=1;
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
function  IndOfMinDistance = FunctionMinDistance(I,J,c1,c2)
%This function calculates the distance between each contour (I,J) point and
%the centre of the contour
%and finds the edge with the  minimum distance 
%      IndOfMinDistance = FunctionMinDistance(I,J,c1,c2)

%      Input Parameters
%      [I,J] points coordinates
%      [c1,c2] coordinates of the centre of the contour

MinDistance = zeros(length(I),1);

MinDistance(:) = sqrt((I(:)-c1).^2+(J(:)-c2).^2);
IndOfMinDistance =find (MinDistance <= min (MinDistance));
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 function [Final_x, Final_y]=PixelsInOrder(x,y,c1,c2);

%Define a cartesian system xy with center[c1,c2],
%Set the points [x,y] in a circular order according to theirs angles 
%       [Final_x, Final_y]=PixelsInOrder(x,y,c1,c2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 c1=round(c1);c2=round(c2);

 vector = [x-c1,y-c2];
 %estimate the angle between the cartesian axis x and the points [x,y]
 Theta = mod(360+180/pi*atan2(vector(:,2),vector(:,1)),360);
 
 %sort the [x,y] points 
 [ThetaSorted,Ind] = sort(Theta);
 
 Final_x=zeros(length(x),1);
 Final_y=zeros(length(y),1);
 Final_x(1:length(x)) =x(Ind);
 Final_y(1:length(y)) =y(Ind);
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function [xi,yi] = Curvefitting(x,y,RES)
% Interpolate(cubic) the points of the image in order to create a close curve
%     [xi,yi] =Curvefitting(x,y,RES)
%
%      RES: resolution desired

%     
%      Modified Code by A. Koulouri
%      
%      Chenyang Xu and Jerry L. Prince, 3/8/96, 6/17/97
%      Copyright (c) 1996-97 by Chenyang Xu and Jerry L. Prince
%      Image Analysis and Communications Lab, Johns Hopkins University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert to column vector
x = x(:); y = y(:);

N = length(x);  

% make it a circular list since we are dealing with closed contour
 x = [x;x(1)];
 y = [y;y(1)];

dx = x([2:N+1])- x(1:N);
dy = y([2:N+1])- y(1:N);
d = sqrt(dx.*dx+dy.*dy);  % compute the distance from previous node for point 2:N+1

d = [0;d];   % point 1 to point 1 is 0 

% now compute the arc length of all the points  from  the first point [x(1),
% y(1)].
% we use matrix multiply to achieve summing 
M = length(d);
d = (d'*uppertri(M,M))';

% now ready to reparametrize the closed curve in terms of arc length
maxd = d(M);

if (maxd/RES<3)
   error('RES too big compare to the length of original curve');
end

di = 0:RES:maxd;

xi = interp1(d,x,di,'cubic');
yi = interp1(d,y,di,'cubic');

N = length(xi);

if (maxd - di(length(di)) <RES/2)  % deal with end boundary condition
   xi = xi(1:N-1);
   yi = yi(1:N-1);
end

function X = uppertri(M,N)
% UPPERTRI   Upper triagonal matrix 
%            UPPER(M,N) is a M-by-N triagonal matrix

[J,I] = meshgrid(1:M,1:N);
X = (J>=I);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%